Einleitung

Dieser Transferbericht untersucht die Einflussfaktoren von tödlichen Unfällen. Er wurde im Rahmen der Vorlesungen “Lineare Regression” und “Zeitreihenanalyse” der BFH im Herbstsemester 2021/22 erfasst, welche von Dr. Raul Gimeno gehalten wurden.

Datenquellen

Unfalldaten

Die Verkehrsunfallstatistik des Kantons Zürich (VUSTA) enthält die Strassenverkehrsunfälle mit Personen- und Sachschäden, die durch die Kantonspolizei Zürich, die Dienstabteilung Verkehr der Stadt Zürich sowie die Stadtpolizei Winterthur registriert wurden. Sie wird einmal jährlich aktualisiert, jeweils gegen Ende des ersten Quartals des Folgejahres.

Zu jedem Strassenverkehrsunfall sind der Unfallort (geokodiert), das Jahr, der Monat, der Wochentag, die Unfallstunde, die Strassenart, der Unfalltyp (ab 1. Juli 2015 inklusive Bagatellunfälle zB. Parkschäden), die Unfallbeteiligung (‘Fussgänger’, ‘Velo (ohne E-Bikes)’, ‘Motorrad’) und die Unfallschwerekategorie verfügbar. Detaillierte Definitionen der Variabeln sind in der Ressource “Minimales Geodatenmodell Strassenverkehrsunfallorte (ASTRA)” dokumentiert.

Die Datei “KTZH_00000718_00001783.csv” wurde am 11.03.2022 heruntergeladen. Am 16.03.2022 wurde eine neue Version veröffentlicht, welche mit dem Script nicht kompatibel ist.

Meteodaten

Der Datensatz umfasst Stundenwerte ab 1992 bis zur letzten aktuellen Stunde aufgeteilt in Jahresdateien. Darin enthalten sind die Stationen Stampfenbachstrasse, Schimmelstrasse und Rosengartenstrasse. Gemessen wird der Luftdruck (p), die Niederschlagsdauer (RainDur), die Globalstrahlung (StrGlo), die Temperatur (T), die relative Luftfeuchtigkeit (Hr), die Windrichtung, die Vektor und Skalar Windgeschwindigkeit. Vor 2018 sind die Skalar Windgeschwindigkeiten aus den 30 Minuten Vektor Daten gerechnet worden. s Die Stundenwerte des laufenden Jahres werden jeweils 30 Minuten nach der vollen Stunde aktualisiert.

Lineare Regression “Unfalldaten”

Daten einlesen

normalizeText <- function(x, length = 9) {
  res <- str_trim(x, side = 'both')
  res <- str_replace_all(res, ' ', '_')
  substr(res, 0, length)
}

# 'basedir' ist eine lokale Variable und verweist auf ein Verzeichnis in meiner Dropbox
accidentsFile <- paste0(basedir, 'KTZH_00000718_00001783.csv')
rawAccident <- read.csv(accidentsFile, header = TRUE, sep = ';', stringsAsFactors = TRUE)

accidentData <- rawAccident %>%
  select(1, AccidentType_de, AccidentSeverityCategory_de, RoadType_de, MunicipalityCode, AccidentYear, AccidentMonth, AccidentWeekDay, AccidentHour) %>%
  mutate(AccidentType_de = normalizeText(AccidentType_de)) %>%
  mutate(AccidentSeverityCategory_de = normalizeText(AccidentSeverityCategory_de, 20)) %>%
  mutate(RoadType_de = normalizeText(RoadType_de)) %>%
  mutate(YearMonth = paste(AccidentYear, sprintf("%02.f", AccidentMonth), sep = '-'))

# Aus Gründen der Übersichtlichkeit wird die UID des Unfalls ausgeblendet
head(accidentData %>% select(-1))
##   AccidentType_de AccidentSeverityCategory_de RoadType_de MunicipalityCode
## 1       Schleuder        Unfall_mit_Sachschad   Hauptstra              247
## 2       Schleuder        Unfall_mit_Sachschad   Nebenstra              261
## 3       Schleuder        Unfall_mit_Leichtver   Nebenstra              261
## 4       Schleuder        Unfall_mit_Sachschad    Autobahn              251
## 5       Schleuder        Unfall_mit_Sachschad      andere              261
## 6       Überquere        Unfall_mit_Leichtver   Nebenstra              261
##   AccidentYear AccidentMonth AccidentWeekDay AccidentHour YearMonth
## 1         2011             1           aw406            0   2011-01
## 2         2011             1           aw406            0   2011-01
## 3         2011             1           aw406            1   2011-01
## 4         2011             1           aw406            1   2011-01
## 5         2011             1           aw406            2   2011-01
## 6         2011             1           aw406            2   2011-01
## Daten einlesen

Unfall nach Typ

accidentByType <- table(accidentData$YearMonth, accidentData$AccidentType) %>%
  as.data.frame() %>%
  pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
  rename(YearMonth = Var1)
tail(accidentByType, 12)
## # A tibble: 12 × 12
##    YearMonth Abbiegeun Andere Auffahrun Einbiegeu Frontalko Fussgänge Parkierun
##    <fct>         <int>  <int>     <int>     <int>     <int>     <int>     <int>
##  1 2020-01          43     18       219        73        13        38       203
##  2 2020-02          64     32       243        79         7        50       205
##  3 2020-03          38     10       147        67         6        23       177
##  4 2020-04          49     12       123        63         7        22       150
##  5 2020-05          70     14       220        86        13        29       189
##  6 2020-06          62     15       249       101        18        48       255
##  7 2020-07          56     22       270        90        23        35       247
##  8 2020-08          56     24       245        73        14        43       257
##  9 2020-09          61     30       261        91        17        43       227
## 10 2020-10          68     23       222       100        22        37       256
## 11 2020-11          69     18       201        70        13        35       204
## 12 2020-12          28     27       195        79        13        38       235
## # … with 4 more variables: Schleuder <int>, Tierunfal <int>, Überholun <int>,
## #   Überquere <int>

Unfall nach Schwerekategorie

accidentBySeverity <- table(accidentData$YearMonth, accidentData$AccidentSeverityCategory) %>%
  as.data.frame() %>%
  pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
  rename(YearMonth = Var1)
tail(accidentBySeverity, 12)
## # A tibble: 12 × 5
##    YearMonth Unfall_mit_Getö… Unfall_mit_Leic… Unfall_mit_Sach… Unfall_mit_Schw…
##    <fct>                <int>            <int>            <int>            <int>
##  1 2020-01                  0              174              946               26
##  2 2020-02                  4              189              977               31
##  3 2020-03                  5              138              740               31
##  4 2020-04                  4              226              629               39
##  5 2020-05                  0              296              910               76
##  6 2020-06                  5              304             1033               66
##  7 2020-07                  2              345             1056               68
##  8 2020-08                  0              299             1048               66
##  9 2020-09                  3              323             1044               71
## 10 2020-10                  2              223             1109               48
## 11 2020-11                  3              198              923               29
## 12 2020-12                  2              166             1009               25

Unfall nach Strassentyp

accidentByRoadType <- table(accidentData$YearMonth, accidentData$RoadType) %>%
  as.data.frame() %>%
  pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
  rename(YearMonth = Var1)
tail(accidentByRoadType, 12)
## # A tibble: 12 × 7
##    YearMonth Autobahn Autostras Hauptstra Nebenanla Nebenstra andere
##    <fct>        <int>     <int>     <int>     <int>     <int>  <int>
##  1 2020-01        173         8       276         2       585    102
##  2 2020-02        182        10       294         1       600    114
##  3 2020-03        121         2       219         0       479     93
##  4 2020-04         72         5       217         3       516     85
##  5 2020-05        132         6       343         2       690    109
##  6 2020-06        150         9       365         2       743    139
##  7 2020-07        185        12       364         4       782    124
##  8 2020-08        191        14       318         1       735    154
##  9 2020-09        197        14       341         4       776    109
## 10 2020-10        162         5       336         5       734    140
## 11 2020-11        148        14       271         2       618    100
## 12 2020-12        142         7       281         0       661    111

Unfall nach Stunde

Wird in der Regressionsanalyse nicht weiter verwendet.

accidentByHour <- table(accidentData$YearMonth, accidentData$AccidentHour) %>%
  as.data.frame() %>%
  pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
  rename(YearMonth = Var1)
tail(accidentByHour, 12)
## # A tibble: 12 × 25
##    YearMonth   `0`   `1`   `2`   `3`   `4`   `5`   `6`   `7`   `8`   `9`  `10`
##    <fct>     <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
##  1 2020-01      24     8     8     6    12    20    33    90    63    54    55
##  2 2020-02      20    15    13    10     3    19    54    77    69    50    80
##  3 2020-03      19    10     5     6     4     5    32    45    64    52    61
##  4 2020-04      15     8     4     4    10     5    13    34    32    38    50
##  5 2020-05      18    11     6     7     9     9    34    48    55    57    73
##  6 2020-06      24    17     4     9     8    12    29    65    70    70    94
##  7 2020-07      20     9     7    14    14    20    45    53    57    78    78
##  8 2020-08      21    17    12    10    12    15    37    60    70    61    96
##  9 2020-09      25    11     9     9    11    20    35    73    76    66    81
## 10 2020-10      16    16    12     4     5    13    50    68    72    65    71
## 11 2020-11      20     8     0     7     5     3    55    51    67    58    64
## 12 2020-12      27    20    11    16     8    22    35    50    63    47    74
## # … with 13 more variables: `11` <int>, `12` <int>, `13` <int>, `14` <int>,
## #   `15` <int>, `16` <int>, `17` <int>, `18` <int>, `19` <int>, `20` <int>,
## #   `21` <int>, `22` <int>, `23` <int>

Unfall nach Wochentag

accidentByWeekDay <- table(accidentData$YearMonth, accidentData$AccidentWeekDay) %>%
  as.data.frame() %>%
  pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
  rename(YearMonth = Var1) %>%
  rename(Mo = aw401, Di = aw402, Mi = aw403, Do = aw404, Fr = aw405, Sa = aw406, So = aw407)
tail(accidentByWeekDay, 12)
## # A tibble: 12 × 8
##    YearMonth    Mo    Di    Mi    Do    Fr    Sa    So
##    <fct>     <int> <int> <int> <int> <int> <int> <int>
##  1 2020-01     170   158   202   202   180   131   103
##  2 2020-02     190   149   177   189   203   205    88
##  3 2020-03     177   141   137   118   141   104    96
##  4 2020-04     117   116   179   158   128   110    90
##  5 2020-05     172   172   195   193   225   189   136
##  6 2020-06     210   280   247   201   192   183    95
##  7 2020-07     193   213   249   256   244   194   122
##  8 2020-08     228   197   198   185   242   196   167
##  9 2020-09     177   252   266   213   229   193   111
## 10 2020-10     166   181   180   236   270   212   137
## 11 2020-11     192   160   188   145   206   152   110
## 12 2020-12     157   235   223   175   190   137    85

Zusammenfassung

Das Datenset erstreckt sich über 120 Monate (Januar 2011 bis Dezember 2022) und enthält rund 150k Unfälle (pro Monat rund 1200).

accidentByMonth <- accidentData %>%
  group_by(YearMonth) %>%
  summarise(TotalAccidents = n())
tail(accidentByMonth, 12)
## # A tibble: 12 × 2
##    YearMonth TotalAccidents
##    <chr>              <int>
##  1 2020-01             1146
##  2 2020-02             1201
##  3 2020-03              914
##  4 2020-04              898
##  5 2020-05             1282
##  6 2020-06             1408
##  7 2020-07             1471
##  8 2020-08             1413
##  9 2020-09             1441
## 10 2020-10             1382
## 11 2020-11             1153
## 12 2020-12             1202
NROW(accidentByMonth)
## [1] 120
mean(accidentByMonth$TotalAccidents)
## [1] 1233.092
sum(accidentByMonth$TotalAccidents)
## [1] 147971
min(accidentByMonth$YearMonth)
## [1] "2011-01"
max(accidentByMonth$YearMonth)
## [1] "2020-12"

Alle Daten zusammenführen

accidentFreq <- accidentByMonth %>%
  # full_join(accidentByHour, by = 'YearMonth') %>%
  full_join(accidentByWeekDay, by = 'YearMonth') %>%
  full_join(accidentByType, by = 'YearMonth') %>%
  full_join(accidentBySeverity, by = 'YearMonth') %>%
  full_join(accidentByRoadType, by = 'YearMonth')

# Sanity checks
head(accidentFreq)
## # A tibble: 6 × 30
##   YearMonth TotalAccidents    Mo    Di    Mi    Do    Fr    Sa    So Abbiegeun
##   <chr>              <int> <int> <int> <int> <int> <int> <int> <int>     <int>
## 1 2011-01             1018   141   112   125   244   160   129   107        41
## 2 2011-02              851   127   110   122   119   145   125   103        34
## 3 2011-03             1071   143   192   161   182   177   127    89        46
## 4 2011-04             1159   149   142   184   195   198   186   105        57
## 5 2011-05             1284   217   203   194   165   190   179   136        70
## 6 2011-06             1127   134   174   213   211   149   152    94        52
## # … with 20 more variables: Andere <int>, Auffahrun <int>, Einbiegeu <int>,
## #   Frontalko <int>, Fussgänge <int>, Parkierun <int>, Schleuder <int>,
## #   Tierunfal <int>, Überholun <int>, Überquere <int>,
## #   Unfall_mit_Getöteten <int>, Unfall_mit_Leichtver <int>,
## #   Unfall_mit_Sachschad <int>, Unfall_mit_Schwerver <int>, Autobahn <int>,
## #   Autostras <int>, Hauptstra <int>, Nebenanla <int>, Nebenstra <int>,
## #   andere <int>
# Still 120 months?
NROW(accidentFreq) == 120
## [1] TRUE

Regressionsanalyse

Hypothese / Untersuchungsgegenstand: Welches sind die Haupteinfluss-Faktoren für “Unfall_mit_Getöteten”?

Mittel: Step-wise Vorwärts- und Rückwärts-Selektion sowie Suche in “beiden Richtungen” mit dem Akaike‘s AIC-Kriterium und dem Schwarz Informationskriterium (BIC).

Vorbereitungen

# Entfernen nicht mehr benötigte Daten
accidents <- accidentFreq %>%
  select(-YearMonth, -TotalAccidents, -Unfall_mit_Sachschad, -Unfall_mit_Leichtver, -Unfall_mit_Schwerver)

# Interzeptmodell als Ausgangmodell
nullModel <- lm(Unfall_mit_Getöteten ~ 1, data = accidents)
fullModel <- lm(Unfall_mit_Getöteten ~ ., data = accidents)

N <- NROW(accidents)

AIC-Kriterium

forwardModelAIC <- step(nullModel,
                        direction = "forward",
                        scope = formula(fullModel), k = 2, trace = 0)

backwardModelAIC <- step(fullModel,
                         direction = "backward",
                         scope = formula(fullModel), k = 2, trace = 0)

bothModelAIC <- step(fullModel,
                     direction = "both",
                     data = accidents, k = 2, trace = 0)

Interpretation

Mit allen drei Verfahren kommen die gleichen Regressoren heraus mit gleichem AIC:

paste(AIC(forwardModelAIC), AIC(backwardModelAIC), AIC(bothModelAIC), sep = ", ")
## [1] "448.006256423627, 448.006256423627, 448.006256423627"
paste(attr(forwardModelAIC$terms, "term"), sep = ", ")
## [1] "Abbiegeun" "Fr"        "Fussgänge" "Tierunfal" "Sa"        "Autobahn" 
## [7] "Frontalko" "So"
paste(attr(backwardModelAIC$terms, "term"), sep = ", ")
## [1] "Fr"        "Sa"        "So"        "Abbiegeun" "Frontalko" "Fussgänge"
## [7] "Tierunfal" "Autobahn"
paste(attr(bothModelAIC$terms, "term"), sep = ", ")
## [1] "Fr"        "Sa"        "So"        "Abbiegeun" "Frontalko" "Fussgänge"
## [7] "Tierunfal" "Autobahn"
finalModelAIC <- lm(Unfall_mit_Getöteten ~
                      Fr +
                        Sa +
                        So +
                        Abbiegeun +
                        Frontalko +
                        Fussgänge +
                        Tierunfal +
                        Autobahn, data = accidents)
summary(finalModelAIC)
## 
## Call:
## lm(formula = Unfall_mit_Getöteten ~ Fr + Sa + So + Abbiegeun + 
##     Frontalko + Fussgänge + Tierunfal + Autobahn, data = accidents)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2681 -0.8868 -0.3055  0.9596  3.5694 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.515703   1.054844   2.385 0.018779 *  
## Fr          -0.015005   0.004392  -3.416 0.000888 ***
## Sa           0.011098   0.005643   1.967 0.051696 .  
## So           0.009357   0.006847   1.367 0.174503    
## Abbiegeun    0.035209   0.013508   2.607 0.010402 *  
## Frontalko    0.057109   0.033439   1.708 0.090458 .  
## Fussgänge    0.035559   0.015551   2.287 0.024112 *  
## Tierunfal   -0.040967   0.021387  -1.916 0.057993 .  
## Autobahn    -0.019077   0.007067  -2.700 0.008028 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.497 on 111 degrees of freedom
## Multiple R-squared:  0.2517, Adjusted R-squared:  0.1977 
## F-statistic: 4.666 on 8 and 111 DF,  p-value: 6.103e-05
layout(matrix(1:4, 2, 2))
plot(finalModelAIC)

jarque.test(finalModelAIC$residuals)
## 
##  Jarque-Bera Normality Test
## 
## data:  finalModelAIC$residuals
## JB = 2.5572, p-value = 0.2784
## alternative hypothesis: greater

Die Residuen sind normalverteilt. Das Modell kann verwendet werden, siehe Analyse / Fazit.

BIC-Kriterium

forwardModelBIC <- step(nullModel,
                        direction = "forward",
                        scope = formula(fullModel), k = log(N), trace = 0)

backwardModelBIC <- step(fullModel,
                         direction = "backward",
                         scope = formula(fullModel), k = log(N), trace = 0)

bothModelBIC <- step(fullModel,
                     direction = "both",
                     data = accidents, k = log(N), trace = 0)

Interpretation

Mit allen drei Verfahren kommen die gleichen Regressoren heraus mit gleichem BIC:

paste(BIC(forwardModelAIC), BIC(backwardModelAIC), BIC(bothModelAIC), sep = ", ")
## [1] "475.881173851447, 475.881173851447, 475.881173851447"
paste(attr(forwardModelAIC$terms, "term"), sep = ", ")
## [1] "Abbiegeun" "Fr"        "Fussgänge" "Tierunfal" "Sa"        "Autobahn" 
## [7] "Frontalko" "So"
paste(attr(backwardModelAIC$terms, "term"), sep = ", ")
## [1] "Fr"        "Sa"        "So"        "Abbiegeun" "Frontalko" "Fussgänge"
## [7] "Tierunfal" "Autobahn"
paste(attr(bothModelAIC$terms, "term"), sep = ", ")
## [1] "Fr"        "Sa"        "So"        "Abbiegeun" "Frontalko" "Fussgänge"
## [7] "Tierunfal" "Autobahn"
finalModelBIC <- lm(Unfall_mit_Getöteten ~ Fr + Abbiegeun, data = accidents)
summary(finalModelBIC)
## 
## Call:
## lm(formula = Unfall_mit_Getöteten ~ Fr + Abbiegeun, data = accidents)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3389 -0.9213 -0.1741  0.7858  4.7462 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.002103   0.773064   2.590 0.010822 *  
## Fr          -0.010850   0.003845  -2.822 0.005611 ** 
## Abbiegeun    0.049589   0.013203   3.756 0.000271 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.584 on 117 degrees of freedom
## Multiple R-squared:  0.1167, Adjusted R-squared:  0.1016 
## F-statistic: 7.732 on 2 and 117 DF,  p-value: 0.0007015
layout(matrix(1:4, 2, 2))
plot(finalModelBIC)

jarque.test(finalModelBIC$residuals)
## 
##  Jarque-Bera Normality Test
## 
## data:  finalModelBIC$residuals
## JB = 6.0343, p-value = 0.04894
## alternative hypothesis: greater

Die Residuen sind nicht normalverteilt. Das Model wird nicht weiter verfolgt.

Analyse und Fazit

Uns ist bewusst, dass mit der vorliegenden Methode keine belastbaren, geschweige denn kausalen Zusammenhänge hergeleitet werden können. Zumal die Kennzahl R2 mit rund 0.25 darauf hinweist, dass das gewählte Regressionsmodell den Regressanden “Unfall_mit_Getöteten” nur ungenügend schwach erklären kann. Dennoch scheint sich zu bewahrheiten, dass die Wochentage Fr (***), Sa (.) und So (nicht signifikant), sowie Abbiegeunfälle (*), Frontalkollisionen (.), Unfälle mit Fussgängern (*) und auf Autobahnen (**) den stärksten Einfluss auf Unfälle mit Todesopfern haben. Dies entspricht unserer Ansicht nach durchaus den Erwartungen. Einzig der Einflussfaktor “Unfälle mit Tieren” (.) kommt etwas unerwartet.

(In Klammern Signifikanzcodes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1)

Das Modell, welches mittels BIC gefunden wurde, wurde verworfen: Die Residuen sind nicht normalverteilt.

Alle Aussagen basieren auf einem Signifikanzniveau von 95%.

Interessant wäre vermutlich die Einflussfaktoren in den Gruppen einzeln zu analysieren, also Wochentage, Unfallart, Strassentyp und ggf. auch Stunde; als auch mit anderen Schwerekategorien zu untersuchen.

Zeitreihenanalyse “Unfall- und Meteodaten”

Daten einlesen

Unfalldaten

Es werden dieselben Daten wie bereits für das Regressionsmodell verwendet.

# function to change column names (remove "Accident", whitespaces and language designator)
process_headers <- function(Headers) {
  NewHeaders <- NULL
  for (i in 1:length(Headers)) {
    tmp <- str_replace(Headers[i], '_', '')
    tmp <- str_replace(tmp, 'Accident', '')
    NewHeaders[i] <- str_replace(tmp, 'en', '')
  }
  return(NewHeaders)
}

# function to extract calendar day from available data
process_day <- function(Month, Weekday) {
  Day <- 1
  for (i in 2:length(Weekday)) {
    # increment day counter if weekday variable changes
    if (Weekday[i] == Weekday[i - 1]) {
      Day[i] <- Day[i - 1]
    } else {
      Day[i] = Day[i - 1] + 1
    }
    # reset if month variable changes
    if (Month[i] != Month[i - 1]) {
      Day[i] = 1
    }
  }
  return(Day)
}

# load raw data into workspace
AccidentFile <- paste0(basedir, 'KTZH_00000718_00001783.csv')
RawAccident <- read.csv(AccidentFile, header = TRUE, sep = ';', stringsAsFactors = TRUE, colClasses = c("AccidentUID" = "character"))

# process raw data to extract information needed for this project
Accident.df <- RawAccident %>%
  select(6, 11, 12, 13, 14, 19, 20, 21, 23, 24, 25, 30, 35) %>%
  drop_na %>%
  rename_with(process_headers) %>%
  mutate(WeekDay = as.numeric(substr(WeekDay, 5, 5))) %>%
  mutate(Day = process_day(Month, WeekDay)) %>%
  mutate(InvolvingPedestrian = ifelse(InvolvingPedestrian == 'true', 1, 0)) %>%
  mutate(InvolvingBicycle = ifelse(InvolvingBicycle == 'true', 1, 0)) %>%
  mutate(InvolvingMotorcycle = ifelse(InvolvingMotorcycle == 'true', 1, 0)) %>%
  mutate(IsSevere = ifelse(grepl('severe', SeverityCategory) | grepl('fatalities', SeverityCategory), 1, 0)) %>%
  mutate(IsDeadly = ifelse(grepl('fatalities', SeverityCategory), 1, 0))

Meteodaten

Die Meteodaten werden direkt von data.stadt-zuerich.ch heruntergeladen und in einen dataframe gespeichert. Wurde die Datei bereits einmal ausgeführt, so ist MeteoDataZurich.Rda bereits verfügbar und wird lokal geladen.

# check if data exists in directory
if (!file.exists(paste0(basedir, 'MeteoDataZurich.Rda'))) {
  # create empty dataframe
  MeteoData.df <- data.frame()
  # get values for 2011 - 2021
  for (i in 2011:2021) {

    url <- paste("https://data.stadt-zuerich.ch/dataset/ugz_meteodaten_stundenmittelwerte/download/ugz_ogd_meteo_h1_", as.character(i), ".csv", sep = '')
    data <- read.csv(url, encoding = 'UTF-8')
    df <- data.frame(data) %>%
      mutate(Year = as.numeric(substr(X.U.FEFF.Datum, 1, 4))) %>%
      mutate(Month = as.numeric(substr(X.U.FEFF.Datum, 6, 7))) %>%
      mutate(Day = as.numeric(substr(X.U.FEFF.Datum, 9, 10))) %>%
      mutate(Hour = as.numeric(substr(X.U.FEFF.Datum, 12, 13))) %>%
      pivot_wider(names_from = Parameter, values_from = Wert) %>%
      select(Year, Month, Day, Hour, T, RainDur) %>%
      group_by(Year, Month, Day, Hour) %>%
      summarise(Temperature = mean(T, na.rm = TRUE), RainDuration = mean(RainDur, na.rm = TRUE)) %>%
      ungroup %>%
      select(Year, Month, Day, Temperature, RainDuration) %>%
      group_by(Year, Month, Day) %>%
      summarise(Temperature = mean(Temperature, na.rm = TRUE), RainDuration = sum(RainDuration, na.rm = TRUE))
    MeteoData.df <- rbind(MeteoData.df, df)
  }
  # save data to Rda file
  save(MeteoData.df, file = paste0(basedir, 'MeteoDataZurich.Rda'))
} else {
  load(paste0(basedir, 'MeteoDataZurich.Rda'))
}

Daten zusammenführen

# add meteo information to dataframe
Accident.df <- merge(Accident.df, MeteoData.df, by = c('Year', 'Month', 'Day'))

Daten bearbeiten und gruppieren

Tagesdaten

# group data into daily sets
AccidentDaily.df <- Accident.df %>%
  select(Year, Month, Day, InvolvingPedestrian, InvolvingBicycle, InvolvingMotorcycle, IsSevere, IsDeadly, Temperature, RainDuration) %>%
  group_by(Year, Month, Day) %>%
  mutate(Total = n()) %>%
  summarise(Total = mean(Total),
            InvolvingPedestrian = sum(InvolvingPedestrian),
            InvolvingBicycle = sum(InvolvingBicycle),
            InvolvingMotorcycle = sum(InvolvingMotorcycle),
            IsSevere = sum(IsSevere),
            IsDeadly = sum(IsDeadly),
            Temperature = mean(Temperature),
            RainDuration = mean(RainDuration)) %>%
  mutate(Date = as.Date(paste(as.character(Year), '-', as.character(Month), '-', as.character(Day), sep = ''), format = '%Y-%m-%d')) %>%
  select(Date, everything()) %>%
  ungroup

Monatsdaten

# further group data into monthly sets
AccidentMonthly.df <- AccidentDaily.df %>%
  group_by(Year, Month) %>%
  summarise(Total = sum(Total),
            InvolvingPedestrian = sum(InvolvingPedestrian),
            InvolvingBicycle = sum(InvolvingBicycle),
            InvolvingMotorcycle = sum(InvolvingMotorcycle),
            IsSevere = sum(IsSevere),
            IsDeadly = sum(IsDeadly),
            Temperature = mean(Temperature),
            RainDuration = sum(RainDuration)) %>%
  mutate(Date = as.Date(paste(as.character(Year), '-', as.character(Month), '-01', sep = ''), format = '%Y-%m-%d')) %>%
  select(Date, everything()) %>%
  ungroup

Grafische Datenexploration

Histogramme

# accidents vs. hour
hist(Accident.df$Hour)

hist(Accident.df$Hour[Accident.df$InvolvingPedestrian == TRUE])

hist(Accident.df$Hour[Accident.df$InvolvingBicycle == TRUE])

hist(Accident.df$Hour[Accident.df$InvolvingMotorcycle == TRUE])

# accidents vs. weekday
hist(Accident.df$WeekDay)

hist(Accident.df$WeekDay[Accident.df$InvolvingPedestrian == TRUE])

hist(Accident.df$WeekDay[Accident.df$InvolvingBicycle == TRUE])

hist(Accident.df$WeekDay[Accident.df$InvolvingMotorcycle == TRUE])

# accidents vs. month
hist(Accident.df$Month)

hist(Accident.df$Month[Accident.df$InvolvingPedestrian == TRUE])

hist(Accident.df$Month[Accident.df$InvolvingBicycle == TRUE])

hist(Accident.df$Month[Accident.df$InvolvingMotorcycle == TRUE])

Korrelationen

# daily data
CorrDataDaily.df <- AccidentDaily.df %>%
  select(Total, InvolvingPedestrian, InvolvingBicycle, InvolvingMotorcycle, IsSevere, IsDeadly, Temperature)

CorrMatrixDaily <- cor(CorrDataDaily.df)

ggcorrplot(CorrMatrixDaily, type = "lower", lab = TRUE)

#ggpairs(CorrDataDaily.df) # attention -> takes a very long time to visualize (large dataset)

# monthly data
CorrDataMonthly.df <- AccidentMonthly.df %>%
  select(Total, InvolvingPedestrian, InvolvingBicycle, InvolvingMotorcycle, IsSevere, IsDeadly, Temperature, RainDuration)

CorrMatrixMonthly <- cor(CorrDataMonthly.df)

ggcorrplot(CorrMatrixMonthly, type = "lower", lab = TRUE)

ggpairs(CorrDataMonthly.df)

Zeitreihen

Zeitreihen plotten

# select some data to make ready to plot with ggplot
AccidentMonthly.plot.df <- AccidentMonthly.df %>%
  select(Date, InvolvingPedestrian, InvolvingBicycle, InvolvingMotorcycle) %>%
  pivot_longer(!Date, names_to = 'Involved', values_to = 'Values')

# plot data using ggplot
ggplot(AccidentMonthly.df, aes(x = Date, y = Total)) +
  geom_line() +
  ylab('Number of Accidents') +
  ggtitle('Total Number of Monthly Accidents') +
  theme_light()

ggplot(AccidentMonthly.plot.df, aes(x = Date, y = Values, color = Involved)) +
  geom_line() +
  ylab('Number of Accidents') +
  ggtitle('Monthly Number of Accidents involving pedestrians, bicycles and motorcycles') +
  theme_light()

Zeitreihenanalyse

# create timeseries from monthly data
InvolvingBicycleMonthly.ts <- ts(AccidentMonthly.df$InvolvingBicycle, start = c(AccidentMonthly.df$Year[1], AccidentMonthly.df$Month[1]), frequency = 12)
InvolvingMotorcycleMonthly.ts <- ts(AccidentMonthly.df$InvolvingMotorcycle, start = c(AccidentMonthly.df$Year[1], AccidentMonthly.df$Month[1]), frequency = 12)

# decompose time series objects
InvolvingBicycleMonthly.ts.m <- decompose(InvolvingBicycleMonthly.ts, type = 'multiplicative')
plot(InvolvingBicycleMonthly.ts.m)

InvolvingMotorcycleMonthly.ts.a <- decompose(InvolvingMotorcycleMonthly.ts, type = 'additive')
plot(InvolvingMotorcycleMonthly.ts.a)

Analyse und Fazit

Grundsätzlich entspricht die Zeitreihenanalyse den Erwartungen, welche man aufgrund der Saisonalität von Unfällen mit Fahrrädern und Motorrädern haben könnte: In den wärmeren Monaten, insbesonderen an warmen Tagen (starke Korrelation zwischen Temperatur und Anzahl Velo- und Motorradunfällen), sind die Unfallzahlen mit Velos und Motorrädern stark erhöht gegenüber den Wintermonaten. Die Zahl der Velounfälle scheint auch mit einem gewissen Trend zu steigen.

Zudem scheint eine negative Korrelation mit Regentagen und Unfällen mit Velos und Motorrädern zu bestehen, was ebenfalls Sinn machen würde.

Im Allgemeinen sieht man nichts Aussergewöhnliches, als das was man aufgrund der Jahreszeiten sowieso vermuten würde.